ESSENTIALLY NON-OSCiLLATORY SHOCK CAPTURING METHODS 
APPLIED TO TURBULENCE AMPLIFICATION IN SHOCK WAVE CALCULATIONS 
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ENO (essentially non-oscillatory) schemes can provide uniformly high order accuracy right 
up to discontinuities while keeping sharp, essentially non-oscillatory shock transitions. Recently 
we obtained an efficient implementation of ENO schemes based on fluxes and TVD Runge-Kutta 
time discretizations. The resulting code is very simple to program for multi-dimensions. ENO 
schemes axe especially suitable for computing problems with BOTH discontinuities AND fine 
structures in smooth regions, such as shock interaction with turbulence, for which results for 
one dimensional and two dimensional Euler equations are presented. We observe much better 
resolution by using third order ENO schemes than by using second order TVD schemes for such 
problems. 

Efficient Implementation of ENO Schemes 


The solutions to systems of hyperbolic conservation laws of the type 

d 

u, + £fi(u) s< =0 (or =gf(u,x,t), a forcing term) (1.1a) 

i=i 

u(x,0) = u°(x) (1.1b) 

where u = ui,...,u m ) T , x = (s 1 , . . . , z rf ), and for real £ = (fi ,...,&*), the combination 
&T& 18 assume( i to have m real eigenvalues and a complete set of eigenvectors, may de- 
velop discontinuities (shocks, contact discontinuities, etc.) regardless of the smoothness of the 
initial condition. Examples of (1.1) include Euler equations of gas dynamics. ENO schemes, 
originally constructed by Haxten, Osher, Engquist and Chakravarthy [1-4], use a local adaptive 

1 Research supported by NSF Grant No. DMS85-03294, DARPA Grant in the ACMP Program, 
ONR Grant N00014-86-K-0691, NASA Langley Grant NAG- 1-270 
2 Research supported by NSF Grant No. DMS88-10150 



stencil to obtain information automatically from regions of smoothness when the solution devel- 
ops discontinuities. As a result, approximations using these methods can obtain uniformly high 
order accuracy right up to discontinuities, while keeping a sharp essentially non-oscillatory shock 
transition. The original ENO schemes in [1-4] used a cell-average framework which involved a 
reconstruction procedure to recover accurate point values from cell averages, and a Lax-Wendroff 
procedure (replacing time derivatives by space derivatives, using the P.D.E.) for the time dis- 
cretization. This can become a bit complicated for multi-dimensional problems [1]. For ease 
of implementation we constructed [7, 8] ENO schemes applying the adaptive stencil idea to the 
numerical fluxes and using a TVD Runge-Kutta type high order time discretization. These ENO 
schemes skip the reconstruction step and the Lax-Wendroif time discretization procedure, hence 
the resulting code is simple for multi-space dimensional problems. 

Let us describe our scheme first in scalar, one dimensional case (d = m = 1 in (1.1)). The 
scheme, in its method-of-lines form, is 




( 1 . 2 ) 


where the numerical flux fj+^ approximates h(xj + ^) to a high order, with h(x) defined by 

/(«(*)) = ^ h(()dS (1.3) 

We first obtain the primitive function of h(x): 

B(x) = f &(0<£e (1.4) 

J — OO 


at x, + j by 



(1.5) 


then construct polynomials interpolating in an ENO fashion, i.e. by obtaining a locally 

“smoothest” stencil starting from one or two points, then adding one point to the stencil at each 
stage by comparing two divided differences and choosing the one which is smaller in absolute 

Ak 

value. fj + 1 is then taken as the derivative of this interpolating polynomial evaluated at x J+ i . 
“Upwinding” is achieved by the initial choice in the stencil-choosing process, and it is also crucial 
for the evident stability of these methods. We also need an entropy fix in any “expansion shock 
cell”. For details, see [7, 8]. 



The time discretization of (1.2) is implemented via a class of TVD Runge-Kutta type methods 
[7]. For example, the third order case is 


: (1) * u (0) + AtL(u (0) ) 

(1.6a) 

+ jAt£(u^) 
4 4 4 

(1.6b) 

m m 1(0) + 2 <2 ) + |a *£(„<’>) 

3 3 3 

(1.6c) 


(1.6d) 


This class of Runge-Kutta methods was shown to have the property that the total variation of 
the spatial part is not increased during the time discretization under a suitable restriction on 

For multi-dimensions the right-hand-side of (1.2) is applied to each of the terms fi(u) Xi in 
(1.1a), keeping all other variables fixed. The Runge-Kutta methods such as (1.6) can still be 
applied. 

For nonlinear systems, we simply apply the algorithms in each local characteristic field. We 
take an 1-dimensional system to exemplify this process. Let A J+ i be some “average” Jacobian 
at Xj + ^ . Examples include Aj + ^ = fj| u _i( Ui+Uj . +1 ) or, in the case of Euler equations of gas 
dynamics Aj + ± = f£| u _ u (*°«) where is the Roe average of u , and u ; +i [6]. We then use 

the eigenvalues of Aj + 1, project to the local characteristic fields: and finally apply our scalar 
algorithms in each of these fields. See [8] for more details. 


2. Numerical Tests — Shock Interaction with Turbulence 

Example 1. We start with one dimensional Euler equations of gas dynamics for a polytropic gas, 
i.e. (1.1) with d = 1, m — 3, and 

u = (p, M, £) t , f(u) = gu + (0, P, qP) T (2.1a) 

where 

P = (7 - 1)(£ - \/>i\ M = pq (2.1b) 

We use 7 » 1.4, and an initial condition 

f p = 3.857143; q = 2.629369; P = 10.333333 when x < -4 
\ p = 1 + £sin5x; q — 0; P = 1 when x > — 4 

If e = 0, this is a pure Mach = 3 shock moving to the right. 


For a detailed linearized analysis see [5]. This linearized analysis, predicts fine structures for 
the density profile because of the different propagation speeds of entropy and acoustic waves. For 


e small (say e — 0.05) we observe results close to linearized analysis. For e = 0.2 we can observe 
nonlinear effects such as additional small shocks in the density profile. 

This is a good test problem because both shocks and fine structures in smooth regions exist. 
Traditional high order methods will develop oscillations near shocks, and TVD methods, while 
nonlinearly stable, will lose resolution for the fine structures because of the degeneracy to first 
order accuracy at smooth critical points. 

In Figure 1-4, the solid lines are numerical solutions of third order ENO scheme (henceforth 
shortened to ENO-3) with 1600 grid points. This can be regarded as a converged solution. From 
Figure 1, we see that ENO-3 with 400 points almost gives a converged solution, while TVD-2 (a 
second order MUSCL type TVD scheme) with 800 points just has roughly the same resolution 
as ENO-3 with 200 points. On the other hand, the improvement of ENO-3 over TVD-2 is not so 
significant for the velocity and pressure profiles (Figure 2), because they both lack any detailed 
structure. 

To further exemplify the advantage of higher order methods, we increase the spatial order of 
our ENO scheme and compare density and entropy profiles with 300 grid points using ENO-3, 4, 
5, 6. We clearly observe better resolution by going to higher spatial orders (Figure 3). In Figure 
3 the time discretization is third order (1.6) with A t decreased for high spatial orders. When 
we use higher order time discretizations as well we observe further improvements in resolution 
(pictures not iiicluded). 

We finally test the effect of physical viscosities by solving the Navier-Stokes equation, i.e. 
(1.1)- (2.1) with a right-hand-side 

4 1 1 2, 2 ' l i I ^ 1 2\ 

’ 3 ' Ife 5 *” rS ' 3 W + ( 7 - lJ-P.-Re-Af'V " 2 9 
We used Pr = 1, M — 3 and gradually increased the Reynolds number Re.- Clearly we observe 
(Figure 4) convergence to Euler’s result as the physical viscosity goes to zero (Re — ► oo). To verify 
the theory (rigorously proven by Kreiss) that for wave lengths > c • the problem is viscosity 
dominated and otherwise essentially inviscid, we re-ran our result with a different frequency for 
the sine wave. We do observe the correctness of the above theory with c«3 in our scaling. The 
pictures axe omitted. 

Example 2. Next we come to two dimensional Euler equations, i.e. (1.1) with d = 2, m = 4, 
and (we use f , g, x, y instead of fi , f 2 , x\ , X 2 ): 

u = (p, M X ,M V , E) T , f(u) = q*u + (0, P, 0, q x P) T 

g(u) = g*u + (0,0,P,g,P) r 



(2.4a) 


where 

P = ( T - 1)(£ - \tf), i 2 =ql + il, M,=pq., M, = pq y (2.4b) 

The test problem we choose is a moving shock interacting with compressible turbulence [9, 
10]. At t = 0, a Mach 8 shock at x = —1.0 is moving into a state with Pr = 1, Pr — 1 
and q x = — §?■ sin 8 r cos (xIcr cos Or + ykR sin Or), q y = cos Or cos (xIcr cos Or + ykR sin Or) 
where Isr = 27r, Or = f, and Cr = We display the results at t — 0.20 in Figure 5. Notice 

that in [9, 10] similar results were obtained using a shock-fitting rather than a shock capturing 
method. This is actually a two dimensional analogue of Example 1 - a combination of shocks and 
fine structures in smooth regions. Hence it is again a good test problem for the high order ENO 
schemes. The successful computation of this example shows that ENO schemes have excellent 
potential for shock-turbulence computations. 
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1(a): ENO-3 with 200 points 
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1(b): ENO-3 with 400 points 
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1(c): TVD-2 with 800 points 


Figure 1: density 
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2(d): TVD-2, pressure 


Figure 2: 200 points 
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5(a): 60 x 40 points, pressure 
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5(d): 120 x 80 points, pressure 



5(b): 60x 40 points, vorticity 


5(e): 120 x 80 points, vorticity 



5(c): 60 x 40 points, entropy 


5(f): 120 x 80 points, entropy 


Figure 5: ENO-3 





